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We perform new longterm (15-16 orbits) simulations of coalescing binary neutron stars in numer¬ 
ical relativity using an updated Einstein’s equation solver, employing low-eccentricity initial data, 
and modeling the neutron stars by a piecewise polytropic equation of state. A convergence study 
shows that our new results converge more rapidly than the third order and using the determined 
convergence order, we construct an extrapolated waveform for which the estimated total phase error 
should be less than 1 radian. We then compare the extrapolated waveforms with those calculated 
by the latest effective-one-body (EOB) formalism in which the so-called tidal deformability, higher 
post-Newtonian corrections, and gravitational self-force effects are taken into account. We show 
that for a binary of compact neutron stars with their radius 11.1km, the waveform by the EOB 
formalism agrees quite well with the numerical waveform so that the total phase error is smaller 
than 1 radian for the total phase of ~ 200 radian up to the merger. By contrast, for a binary of less 
compact neutron stars with their radius 13.6 km, the EOB and numerical waveforms disagree with 
each other in the last few wave cycles, resulting in the total phase error of ~ 3 radian. 

PACS numbers: 04.25.D-, 04.30.-w, 04.40.Dg 


I. INTRODUCTION 

The inspiral and merger of coalescing compact bina¬ 
ries are among the most promising sources for ground- 
based kilometer-size laser-interferometric gravitational- 
wave detectors [l|. A statistical study based on the 
stellar-evolution synthesis (e.g., Refs. @,131) suggests that 
the detection rate for them will be 1-100 yr“^ for 
advanced detectors, i.e., advanced LIGO [i|, advanced 
VIRGO Q , and KAGRA @ , which will sequentially start 
operation in the coming years. 

One of the important steps after the first detection of 
gravitational waves from binary neutron stars (and also 
a black hole-neutron star binary) will be to extract bi¬ 
nary parameters such as mass, spin, and radius of each 
object in the binary systems. In particular, the mass and 
radius (or a quantity related to it) of the neutron stars 
have invaluable information for determining the equation 
of state (EOS) of the neutron-star matter, which is still 
poorly known. The mass of two neutron stars will be de¬ 
termined with a high accuracy < 1%, if the gravitational- 
wave signals in the inspiral stage are detected with the 
signal-to-noise ratio ^10 and the neutron-star spins are 
supposed to be negligible Q- On the other hand, de¬ 
termining the parameters related to the neutron-star ra¬ 
dius is the challenging issue although it has to be done 
for constraining the neutron-star EOS (e.g.. Ref. IM). 
Among other possible methods, extracting the tidal de¬ 
formability of the neutron stars from gravitational waves 
emitted binary-neutron-star inspirals is one of the most 
promising methods [HHil- For employing this method, 


we have to prepare a theoretical template of gravitational 
waves from binary-neutron-star inspirals taking into ac¬ 
count tidal-deformation effects that influence the dynam¬ 
ics of the late inspiral orbits (e.g.. Ref. H)- Hence, the¬ 
oretically deriving a precise gravitational waveform for 
binary-neutron-star inspirals including the tidal effects is 
an urgent task. 

A post-Newtonian (PN) gravitational waveform for the 
early stage of binary-neutron-star inspirals (with the fre¬ 
quency / < 400 Hz) was first derived by Flanagan and 
Hinderer including the leading-order tidal effects [l^ . 
They showed that the tidal effect for the evolution of the 
gravitational-wave phase could be described only by the 
tidal deformability of neutron stars. They also found that 
the tidal deformability of neutron stars could be mea¬ 
sured by the advanced gravitational-wave detectors by 
using the gravitational-wave signals for / = 10 - 400 Hz, 
if the tidal deformability of neutron stars is sufficiently 
large or if we could observe an event with a high signal- 
to-noise ratio (see also Ref. [l^). If the waveform is ex¬ 
tended to the higher frequency range, the measurability 
can be significantly improved. In the PN approach, how¬ 
ever, the uncertainty of the higher-order PN terms pre¬ 
vents us to construct the accurate waveform at the higher 
frequency [l6l - [l8j| . 

To overcome the ambiguity in the higher PN terms, an 
effective one body (EOB) formalism with the tidal effects 
has been explored [11, [H, HI . In this approach, the non- 
tidal part is calibrated using the results of binary-black- 
hole simulations. Damour and his collaborators [2l| 
subsequently explored the measurability of the tidal de- 
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formability with the advanced gravitational-wave detec¬ 
tors employing an EOB formalism including tidal effects 
up to the second PN order. They concluded that the 
tidal deformability of neutron stars could be measured by 
the advanced gravitational-wave detectors if the signal- 
to-noise ratio of the gravitational-wave signal is higher 
than ^ 16 for any EOS that satisfies the constraint of 
the maximum gravitational mass > 2Mq The key 

assumption of their study is that the EOB approach is 
valid up to the onset of the merger of binary neutron 
stars. However, in the stage just before the merger, ef¬ 
fects such as nonlinear tidal-deformation effects, which 
are not taken into account in the current EOB formal¬ 
ism, could come into play (see, e.g., Ref. M)- 

For precisely understanding the orbital motion and 
the waveform in the late inspiral stage of binary neu¬ 
tron stars, a high-resolution numerical-relativity (NR) 
simulation with appropriately physical setting is obvi¬ 
ously necessary. Recently, long-term simulations for 
binary-neutron-star inspirals were performed by three 
groups [iO, [23 - [^ aiming at the derivation of accurate 
gravitational waveforms for the late inspiral stage. They 
followed the late binary inspiral for < 10 orbits up to the 
onset of the merger. However, in their numerical simu¬ 
lations, an unphysical residual eccentricity is present in 
the initial data. This seriously made their results less 
accurate, because binary neutron stars in the late inspi¬ 
ral stage are believed to have a quasi-circular orbit with 
negligible eccentricity. In the present work, we simu¬ 
late binary-neutron-star inspirals for a longer term with 
more physical initial data in which the eccentricity is suf¬ 
ficiently small (less than 10“^) In addition, we per¬ 
form the simulations employing a formalism in which the 
constraint violation can be suppressed to a level much 
smaller than that in our previous study [^ . As a result, 
we can obtain an extrapolated waveform in a much more 
accurate and reliable manner than in our previous study. 

The paper is organized as follows. In Sec. H, we sum¬ 
marize the formulation and numerical schemes employed 
in our numerical-relativity study, and also review the 
EOS employed. In Sec. HI, we describe our method for 
deriving an extrapolated gravitational waveform, show¬ 
ing the resulting waveforms that are much more accurate 
than those derived in our previous study [2^. We then 
compare our extrapolated waveforms with those derived 
by the latest EOB approach and examine its accuracy in 
Sec. IV. Section V is devoted to a summary. Throughout 
this paper, we employ the geometrical units of c = G = 1 
where c and G are the speed of light and the gravitational 
constant, respectively. 


^ We note that R. Haas and his collaborators (SXS collaboration) 
have also derived the waveforms of small eccentricity in their 
longterm simulations, although their results have not been pub¬ 
lished yet. 


II. FORMULATION FOR 
NUMERICAL-RELATIVITY SIMULATION 

In this section, we briefly describe the formulation and 
the numerical schemes of our numerical-relativity simu¬ 
lation employed in this work. 

A. Evolution and Initial Condition 

We follow the inspiral and early stage of the merger 
of binary neutron stars using our numerical-relativity 

code^ _ SACRA, for which the details are described in 

Ref. [291 . In this work, we employ a moving puncture ver¬ 
sion of the Baumgarte-Shapiro-Shibata-Nakamura for¬ 
malism [30l| . locally incorporating a Z4c-type constraint 
propagation prescription |3ll | (see for our implemen¬ 
tation) for a solution of Einstein’s equation. The con¬ 
straint propagation from the neutron-star’s outer region 
plays a crucial role for reducing the constraint viola¬ 
tion and for improv ing the order of the convergence as 
discovered in Ref. |3l| . In our numerical simulation, 
a fourth-order finite differencing scheme in space and 
time is used implementing an adaptive mesh refinement 
(AMR) algorithm. At refinement boundaries, a second- 
order interpolation scheme is partly used. The advection 
terms are evaluated b y fo urth-order lop-sided upwind- 
type finite differencing . A fourth-order Runge-Kutta 
method is employed for the time evolution. For the hy¬ 
drodynamics, a high-resolution central scheme based on 
a Kurganov-Tadmor scheme with a third-order piece- 

wise parabolic interpolation and with a steep min-mod 
limiter is employed. 

In this work, we prepare nine refinement levels for the 
AMR computational domain. Specifically, two sets of 
four finer domains comoving with each neutron star cover 
the region of their vicinity. The other five coarser do¬ 
mains cover both neutron stars by a wider domain with 
their origins fixed at the center of the mass of the binary 
system. Each refinement domain consists of a uniform, 
vertex-centered Cartesian grid with (2V-|-1, 2N -|-1, V-|- 
1) grid points for (x, y, z) (the equatorial plane symme¬ 
try at z = 0 is imposed). The half of the edge length 
of the largest domain (i.e., the distance from the origin 
to outer boundaries along each axis) is denoted by L, 
which is chosen to be larger than Aq, where Aq = tt/Oq 
is the initial wavelength of gravitational waves and Hq is 
the initial orbital angular velocity. The grid spacing for 
each domain is Axi = L/{2^N), where Z = 0 — 8. In this 
work, we choose N = 72, 60, 48, and 40 for examining 
the convergence properties of numerical results. With 
the highest grid resolution (for N = 72), the semimajor 
diameter of each neutron star is covered by about 120 
grid points. 

We prepare binary neutron stars in quasi-circular or¬ 
bits for the initial condition of numerical simulations. 
These initial conditions are numerically obtained by us¬ 
ing a spectral-method library, LORENE [35|. We fol- 
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TABLE I: Equations of state (EOS) employed, the radius and the tidal Love number of Z = (2, 3, 4) of spherical neutron stars 
of mass 1.35Mq, the radius of light ring orbit, angular velocity of initial data, and the finest grid spacing in the four different 
resolution runs, mo denotes the total mass of the system. In this study, it is 2.7Mq. 


EOS 

R 1.35 (km) 

ZC2,1.35 

^3,1.35 

ZC4,1.35 

rLR 

mofZo 

A^s (km) 

APR4 

11.1 

0.0908 

0.0234 

0.00884 

3.61 

0.0156 

0.140, 0.167, 0.209, 0.251 

H4 

13.6 

0.115 

0.0326 

0.0133 

4.21 

0.0155 

0.183, 0.220, 0.274, 0.329 


low 15-16 orbits in this study. To do so, the orbital 
angular velocity of the initial configuration is chosen 
to be mof^o ~ 0.0155 (/ « 370 Hz for the total mass 
mo = 2.7Mq, i.e., each mass of neutron stars is 1 . 35 M 0 ). 
The neutron stars are assumed to have an irrotational 
velocity field, which is believed to be an astrophysically 
realistic configuration [H, [s^. The parameters for the 
initial models are listed in Table ID 

For the computation of an accurate gravitational wave¬ 
form in numerical simulations, we have to employ initial 
data of a quasi-circular orbit of negligible eccentricity. 
Namely, the eccentricity of the initial binary orbit has to 
be reduced to be as small as possible. Such initial data 
are constructed by an eccentricity-reduction procedure 
described in [s^- For the initial data employed in this 
work, the residual eccentricity is < 10 “^. 


parts as 

H = Pcoid(p) + Tth , £ = e^coid(p) + £th- (2.2) 

The cold parts of both variables are calculated using the 
original piecewise polytropic EOS from p, and then the 
thermal part of the specific internal energy is defined 
from £ as eth = e — ecoid(p)- Because Sth vanishes in 
the absence of shock heating, it is regarded as the finite- 
temperature part determined by the shock heating in the 
present context. For the thermal pressure, a F-law ideal- 
gas EOS was adopted as 

Pth = (Fth - l)peth- (2.3) 

Eollowing our latest works [H, |4l|, Fth is chosen to be 

1 . 8 . 


B. Equation of State 

Following previous works [l^, |32| , we employ a param¬ 
eterized piecewise-polytropic equation of state proposed 
by Read and her collaborators [^. This EOS is written 
in terms of four segments of polytropes 

P = Kip^' ( for Pi < p < pi+i, 0 < i < 3), (2.1) 

where p is the rest-mass density, P is the pressure, 
Ki is a polytropic constant, and F^ is an adiabatic in¬ 
dex. At each boundary of the piecewise polytropes, 
p = Pi, the pressure is required to be continuous, i.e.. 
Following Read and her collabo¬ 
rators, these parameters are determined in the follow¬ 
ing manner [^: The crust EOS is fixed by setting 
Fq = 1.3562395 and Kq = 3.594 x 10^^ in cgs units. The 
values of the boundary density is set as p 2 = g/cvo? 

and p 3 = 10^® °g/cm^. With this preparation, the 
following four parameters become free parameters that 
should be given; {Pi,Fi,F2,F3}. Here, Pi is the pres¬ 
sure at p = p 2 , and for a given value of this, Ki and Ki 
are determined by Ki = Pi/p 2 ^ and ATi+i = ■ 

In this work, we choose two sets of piecewise-polytropic 
EOS mimicking APR4 and H4 |40l| EOS (see Table 
1 of Ref. [dlj for the four parameters). 

In numerical simulations, we employ a modified version 
of the piecewise polytropic EOS to approximately take 
into account thermal effects, which play a role in the 
merger phase. In this EOS, we decompose the pressure 
and specific internal energy into the cold and thermal 


C. Extraction of Gravitational waves 


Gravitational waves are extracted from the outgoing- 
component of complex Weyl scalar ^'4 [ 2 ^. From this, 
gravitational waveforms are determined in spherical co¬ 
ordinates (r, P, (j)) by 


h := h_|_(Z) — ihx{t) = — lim [ dt' [ dt"'i> 4 {t", r). 


(2.4) 


Here, we omit arguments 6 and </>. 414 can be expanded 
in the form 

4 / 4 (Z, r, eA)=J2 r)_2Yim{0. 4>), (2.5) 

Im 

where - 2 Yim denotes the spin-weighted spherical har¬ 
monics of weight —2 and are expansion coefficients 
defined by this equation. In this work, we focus only 
on the (Z, |m|) = (2,2) mode because we pay attention 
only to the equal-mass binary, and hence, this quadrupole 
mode is the dominant one. 

We evaluate 4'4 at a finite spherical-coordinate radius 
r/mo = 100-240. To compare the waveforms extracted 
at different radii, we use the retarded time defined by 

( 2 . 6 ) 


where r* is the so-called tortoise coordinate defined by 

n := TA-I- 2moln - 1^ , (2.7) 

with ta := A/Att and A the proper area of the extrac¬ 

tion sphere. 
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III. RECIPE FOR CONSTRUCTING AN 
EXTRAPOLATED WAVEFORM 

In this section, we present our prescription for deriving 
an extrapolated gravitational waveform from raw numer¬ 
ical data of 4 / 4 , and show that the resulting waveforms 
have a good accuracy that can be compared carefully 
with the EOB results. 


A. Extrapolation to infinite extraction radius 


As we mentioned in the previous section, we extract 'I '4 
at several coordinate radii, 100-240mo, and then, this 
complex Weyl scalar is decomposed into the spherical 
harmonics components, 'P 4 ™. Since the waveform of 'I '4 
extracted at a finite radius, ro, is systematically different 
from that at null infinity, we first compute an extrapo¬ 
lated waveform at tq —^ 00 using the Nakano’s method 

as [mil 


^'4’"’°°(iretUo) 


= C'(ro) 4 ' 4 '"(tretUo) 


2r-A 


P^ret 


4'4’™(t',ro)dt' , 
(3.1) 


where C'('ro) is a function of vq. Since our coordinates 
are similar to isotropic coordinates of non-rotating black 
holes, we choose ta = r'o[l -|- mo/( 2 ro)]^. C'(r-o) depends 
on the choice of the tetrad components; for our choice, 
it is appropriate to choose C'(ro) = 1 — 2mo/?’A. In this 
setting, tret at r = ro is given by Eqs. ( 1 ^ and (I 2 J 1 ) . 

The left panel of Fig. [T] plots the real part of 
'I' 4 ’^’°°(tretUo) for Several choices of tq. The right 
panel shows the evolution of the absolute amplitude of 
ro). These show that the extrapolated wave¬ 
forms depend very weakly on the extraction radius, ro 
(see Ref. [i^ for the reason). 

We then have to calibrate how weakly the resulting 
extrapolated waveforms, 4 ' 4 ’^’°“(tretUo); depend on ro 
and have to estimate the systematic error in this quan¬ 
tity. We find that the systematic error in phase decreases 
approximately in proportional to (cf. the left lower 
panel of Fig. [T] that indeed shows this property). Fig¬ 
ure [T] implies that for ro > 200 mo, the systematic error 
in phase is smaller than 0.3 radian. This value is smaller 
than the error in the extrapolated waveform finally ob¬ 
tained (associated with the uncertainty in the resolution 
extrapolation), and can be accepted in the present nu¬ 
merical study. Note that this phase error is systematic 
and could be subtracted, although we do not do so in 
this work. 

By contrast, the systematic error in amplitude is ap¬ 
preciable, i.e., 1-2 percents even for ro ^ 200mo. For 
suppressing this error, we might have to enlarge the com¬ 
putational domain for the wave extraction. However, 
this error size is smaller than another error associated 


with the spurious short-term modulation in the numeri¬ 
cal gravitational-wave amplitude as reported in Ref. : 
The right panel of Fig. [T] shows that a modulation in 
the amplitude is present with its fluctuation amplitude 
of < 3% in particular in the early stages of the nu¬ 
merical waveform. Since this error was not able to be 
cleaned up, we do not take a further extrapolation of 
|4'4™’°°(ti.etUo)| for ro — 00 . Thus, in this work, we em¬ 
ploy ro) computed from the data extracted 

at ro = 200too [hereafter written as 'I' 4 ™’°°(ti.et)] with¬ 
out further processing and perform subsequent analyses 
keeping in mind that in the amplitude extrapolated by 
Eq. (13.11) . there could exist a local error in magnitude up 
to ~ 3% of the exact amplitude (note that in average the 
error would be much smaller than 3%). 


B. Extrapolation for zero-grid spacing limit 

Next, we consider the resolution extrapolation for the 
limit Axs —t 0. For this task, numerical simulations have 
to be performed for more than three grid resolutions. In 
this study, we performed four simulations for each model 
employing four different grid resolutions (cf. Table U for 
the finest grid spacing, Accg, for each run). For each run, 
we extracted the numerical waveform at tq = 200too and 
then performed the extrapolation of ro —>■ 00 as described 
in Eq. (13.11) . 

We then need to perform an extrapolation procedure 
of taking the zero grid-spacing limit for obtaining an ap¬ 
proximately exact solution. For this procedure, we first 
analyze the relation of the time to the merger, tmrg, as 
a function of Axs following Ref. [l^. Here, the merger 
time, tmrg, is defined as the time at which the maximum 
value of |'I '4 ^’““(fret)! is recorded. Then, it is found that 
tmrg converges to an unknown exact value at 4th or¬ 
der (see below for more detailed analysis), fmrg is larger 
for the better grid resolutions because for the lower grid 
resolutions, the numerical dissipation is larger and the 
inspiraling process is spuriously accelerated. This numer¬ 
ical error is universally present for finite values of Axg; 
namely, for any inspiraling stage in any numerical sim¬ 
ulations, the error is always present. For obtaining the 
“exact” waveform, thus, we always need an extrapolation 
procedure. Then, the next question is how to extrapo¬ 
late the waveform for the limit Axs —t 0. We propose 
the following method in this study. 

We first determine the gravitational waveform and 
time evolution of the angular frequency as functions of 
^ret by integrating 4 ' 4 ’”’°“(tj.et) for each raw numerical 
data. Here, the gravitational waveform for each multi¬ 
pole mode satisfies [see Eq. (EH)] 

hi,m jg obtained by the double time integration of 

For this procedure, we employ the method of Ref. [dj, 





5 



10 20 30 40 50 

tret (ms) 


O 

s 

o 

Lh 


(N 



0.0004 

0.00035 

0.0003 

0.00025 

0.0002 

0.00015 

0.0001 



Let (ms) 


FIG. 1: The waveform (real part; left) and amplitude (right) of '^^’^’““(ro, fret) as functions of fret for several values of ro for 
the run with H4 EOS and the best grid resolution {N = 72). The lower plot of the left panel shows the phase differences of 
tj>2,2,oo(^o) relative to = 237mo). 


written as 

/ ^.l,m^oo / \ 

doj - - -exp(za;tret), (3.3) 

max(w, Wcut) 

where is the Fourier transform of 

and Wcut is chosen to be I.GOq. (Note that at the initial 
stage, the value of uj is 2r2o > Wcut). We recall again that 
in this paper we pay attention only to / = |m| = 2 modes 
because these are the dominant modes in particular for 
the equal-mass binaries. Then, from Eq. (lO) . we deter¬ 
mine the evolution of the amplitude, i.e., 
as a function of tret- 

Using Eq. (13.31) . we can also define the evolution of the 
angular frequency as 


w(tret) 1^2,2| ’ 

and then, the evolution of the gravitational-wave phase 
is calculated by 


4* (tret) 



dt' u;(t'). 


(3.5) 


phase evolution is spuriously faster for the poorer grid 
resolutions. However, we already know that the merger 
time converges approximately at 4th-order. Taking into 
account this fact, we stretch the time axis for the grav¬ 
itational waveform by an appropriate factor as t —^ r^t 
where t]{> 1) is the constant stretching factor. This fac¬ 
tor should be larger for the results of the poorer grid reso¬ 
lutions. Here, this reprocessing is performed in the same 
manner as in [^ : Uet and <1> are modified as tret ^ viret 
and $ —7- 77$. We will show that the phase evolution 
matches very well among the waveforms with different 
grid resolutions after this scaling performed in terms of 
this single parameter ij. Later, ij will be also used for 
determining the convergence order and for obtaining the 
resolution-extrapolated waveform. 

As a first step for this stretching procedure, we have 
to determine the values of rj. As the first substep, we 
carry out a procedure for finding the minimum value of 
the following integral 


Af 


I = min 
v',4> . 


dtr 


exp [ir]'^ 2 Wtr:et) + *<(>] 


- A^’^(Uet)exp[7$i(Uet)] , (3.7) 


Now, using and $, the quadrupole gravitational 

waveform can be written by 

(l^’^(Uet) = A^’^(Uet) exp [l4>(tret)] • (3.6) 

Eigure [2] plots the resulting gravitational waveforms 
and the evolution of $ obtained in the simulations with 
different grid resolutions for the models with H4 (left) 
and APR4 EOS (right). The upper panels plot the grav¬ 
itational waveforms and these show that the merger time 
is earlier for the poorer grid resolutions. The middle pan¬ 
els plot the integrated wave phases for the pure numer¬ 
ical results with no reprocessing. These show that the 


where A^'^ and are, respectively, the amplitude and 
integrated phase of the gravitational waveform for the 
best-resolved run (N = 72) and A^’^ and $2 are those 
for less-resolved runs. The free parameters, rj' and (j), are 
varied for a wide range and from 0 to 27r, respectively, 
to search for the possible minimum value of I. ti and tf 
are chosen to be 5 ms and tmrg of the best resolved run, 
respectively. Here, the reason for choosing ti = 5 ms is 
that for their early stage with fret 5 ms, the numerical 
waveforms have a relatively large modulation in ampli¬ 
tude and phase due to junk radiation. 

We find for our present simulation results that for the 
second-finest, third-finest, and poorest resolution runs. 
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FIG. 2: The gravitational waveforms and the evolntion of the gravitational-wave phase for four different grid-resolution runs 
with the H4 EOS (left three panels) and with the APR4 EOS (right three panels). N indicates the grid resolution, Aa:8 oc N~^. 
The upper panels show the gravitational waveforms for three different grid resolutions. The middle panels show the pure 
numerical wave phases and the bottom panels show the results obtained after the stretching of time and phase according to 
the convergence property (for N = 40, 48, and 60). The lower plots in middle and bottom panels show the phase disagreement 
between the purely numerical wave phase for N = 12 and the lower-resolution results. Note that for N = 72, tmrg = 58.43 (ms) 
for the H4 EOS and tmrg = 61.08 (ms) for the APR4 EOS. 


T]' = 1.00646, 1.02241, and 1.06000 for the H4 EOS and 
r]' = 1.00650, 1.02931, and 1.09118 for the APR4 EOS. 
The mismatched factors, respectively, are ///q = 7.4 x 
10"®, 2.3 X 10-^ and 1.4 x IQ-^ for the H4 EOS and 
I/h = 7.4 X 10"®, 1.1 X 10-4, and 1.4 x IQ-^ for the 


APR4 EOS. Here, we define 



The cross correlation of two waveforms is approximately 
estimated as 1 — yjl/2 Iq. This implies that the cross 
correlation between the waveforms of the best-resolved 
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run and reprocessed less-resolved runs are approximately 
99.9%, 99.8%, 99.4% for the H4 EOS and 99.9%, 99.5%, 
and 98.2%, respectively. This shows that the accuracy is 
not very good in the low-resolution runs for the APR4 
EOS, for which the compactness is larger than that for 
the H4 EOS, and hence, a finer grid resolution would be 
necessary for a well-resolved simulation. Eor both EOS, 
the reprocessed waveforms in the poorest-resolution run 
are found to be not very accurate, and hence, in the 
following, we will perform a convergence study employing 
the waveforms of the first-, second-, third-resolved runs 
(labeled by = 72, 60, and 48, respectively). 

The bottom panels of Eig. [5]show the results obtained 
for this time-stretching procedure. It is found that four 
curves of <i> originally with different grid resolutions ap¬ 
proximately overlap with each other. In particular, the 
degree of the overlapping is quite good between the finest 
and second-finest runs (see the difference of the inte¬ 
grated phase shown in the lower plot of the bottom pan¬ 
els of Eig. El): For both EOS, the disagreement of $ for 
these reprocessed data is much smaller than 0.1 radian 
except for the final moment of the last orbits, at which 
the disagreement steeply increases: however it is at most 
~ 0.2 radian. This suggests that the time stretching 
method can be used for obtaining the extrapolated wave¬ 
form for Axs —?► 0 if we accept the error of the integrated 
phase up to ~ 0.2 radian. 

We next try to obtain an extrapolated waveform for 
Axg —>■ 0 by using the time stretching method for the 
well resolved models. For this procedure, we have to de¬ 
termine the order of the convergence appropriately. In 
the above, we found that the numerical waveform in the 
poorest run is not very reliable even after the reprocess¬ 
ing. Thus, we determine the order of the convergence 
from the three better-resolved runs. (Note that if we em¬ 
ploy the poorest-resolved waveforms for determining it, 
the order of the convergence is spuriously overestimated.) 
Using the values of rj' — 1, the order of the convergence, 
p, is determined from 


(72/48)P - 1 
(72/60)P- 1 


0.02241 

0.00646 


for H4, 


0.02931 

0.00650 


for APR4, 


(3.9) 


which give p ~ 3.42 and 5.10, for the H4 and APR4 EOS, 
respectively. This indicates that the stretching factor for 
the best-resolved run to reproduce the limiting waveform 
with Axs —0 is r; ss 1.00746 and 1.00424 for the H4 
and APR4 EOS, respectively. This implies that for these 
models, the exact merger time would be tmrg ~ 58.87ms 
and 61.34 ms, respectively, whereas they were 58.43 ms 
and 61.08 ms for the best-resolved run. Namely, the error 
in the merger time is still much larger than 0.1ms even 
for the best-resolved run: For obtaining the waveforms of 
the error in the merger time smaller than 0.1 ms, a simu¬ 
lation with N > 100 would be necessary. We note that if 
we extrapolate the value of rj for the best-resolved runs 


assuming the third- and fourth-order convergences of rj', 
the value of rj becomes, respectively, 1.00887 and 1.00601 
for the H4 EOS and 1.00893 and 1.00605 for the APR4 
EOS. For the hypothetical fourth-order convergence, the 
predicted merger time would be tmrg = 58.78 ms for the 
H4 EOS and 61.45ms for the APR4 EOS. Thus, it is 
safe to keep in mind that the extrapolated merger time 
still has an error of ^ 0.1 ms due to the uncertainty in p. 
Since the merger time is ^ 60 ms and total gravitational- 
wave phase is ^ 200 radian for both EOS, we should keep 
in mind the phase error of 200 x (0.1/60) ^ 0.3 radian. 


IV. COMPARISON BETWEEN 
NUMERICAL-RELATIVITY AND 
EFFECTIVE-ONE-BODY WAVEFORMS 


Figure |3] plots the extrapolated gravitational wave¬ 
forms, the associated frequency, and the integrated 
gravitational-wave phase. For comparison, we plot the 
results by an effective-one-body (EOB) approach Bin, 
Idsjl (see appendix A for the EOB formalism that we em¬ 
ploy in this work). To align the time and phase of the 
numerical and EOB waveforms, we first calculate a cor¬ 
relation like Eq. (1X71) for 5 ms < Uet < 20 ms between 
the numerical and EOB waveforms, varying the time and 
phase of the EOB waveform. These parameters are de¬ 
termined by searching for the set of the values that give 
the minimum of this integral. 

Figure in shows that up to / ~ 700 Hz (at Uet ~ 54 ms), 
the EOB result well reproduces the extrapolated wave¬ 
forms for both H4 and APR4 EOS: In particular for the 
APR4 EOS for which the compactness is large and the 
tidal deformability is small, the agreement is quite good. 
For both EOS, the error in the frequency is smaller than 
1% and the phase error is smaller than 0.1 radian for 
/ 700 Hz (with Uet > 5 ms). However, for the last a few 

cycles, the agreement between the extrapolated and EOB 
waveforms becomes poor. Here, note that this disagree¬ 
ment cannot be explained by the error in the numerical 
waveform, because we have already estimated that the 
phase error in the numerical waveform would be smaller 
than ~ 0.3 radian. The magnitude of the error is larger 
for the H4 EOS. The possible reason for this disagree¬ 
ment is that in the current version of the EOB formal¬ 
ism, the tidal effects are not fully taken into account (e.g., 
non-linear tidal effects and non-stationary effects are not 
included). Namely, if the degree of the tidal deformation 
becomes high, the approximation could be poor. 

In the final inspiraling stage for the model with the H4 
EOS, the neutron stars are significantly deformed, and 
the attractive force associated with the tidal deforma¬ 
tion is enhanced: The relative fraction of the approach¬ 
ing velocity induced by the tidal effect to that by other 
general relativistic effects such as gravitational-radiation 
reaction is larger for the binary of larger-radius neutron 
stars. The missing tidal effects could give a significant 
damage in the current version of the EOB formalism. By 
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FIG. 3: The extrapolated gravitational waveform and related quantities for the models with the H4 (left) and APR4 EOS 
(right). Top: The extrapolated waveforms for the best-resolved {N = 72) and second-best-resolved {N = 60) runs are plotted 
(two waveforms overlap quite well with each other and we cannot distinguish them in the figure). The waveform by an effective- 
one-body calculation is plotted together. The lower panels focus on the late inspiral waveforms. Middle: The extrapolated 
gravitational-wave frequency. In the lower panel of this, the absolute difference between the extrapolated result (with N = 72) 
and EOB result is shown. Bottom: The extrapolated gravitational-wave phase. In the lower panel of this, the difference 
between the extrapolated result and EOB result is shown. We aligned the phases of the extrapolated and EOB waveforms at 
fret = 5 ms. 


contrast, for the model with the APR4 EOS, the agree¬ 
ment between the extrapolated and EOB waveforms is 
quite good even at the last orbit. The total phase error 
is smaller than 0.7 radian, which is comparable to that 
in the error associated with the uncertainty of the extrap¬ 


olation. This implies that for the binary of small-radius 
neutron stars, the current version of the EOB formalism 
would be already robust if we accept the phase error of 
~ 1 radian (see also Ref. [^). 

The missing tidal effects in the EOB formalism cannot 
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tiet (ms) 


FIG. 4: The evolution of the mismatch, 7m(tret), between 
the extrapolated waveform and EOB waveforms. “EOB ±3” 
denotes that the EOB waveforms are employed artificially in¬ 
creasing or decreasing the neutron-star compactness by 3%. 


be compensated even if we artificially modify the value of 
the tidal deformability (or compactness). Figure |4] plots 
the evolution for the degree of the mismatch between 
the extrapolated waveform and the EOB waveform. For 
comparison, we calculated the mismatch employing the 
EOB waveforms in which the compactness of neutron 
stars is varied by ±3%. Here, the mismatch is defined by 


_1 {h /leobl^ ^eob) 

2 (^eobl^eob) 


where 


{hi\h2) 



^l(^ret)^2(^ret)'^^ret 


(4.1) 


(4.2) 


Again, ti is chosen to be 5 ms. Here, h and /leob denote an 
extrapolated waveform and a waveform by the EOB for¬ 
malism, respectively. We note that the following relation 
is approximately satisfied for small values of Im{tret)' 


1 .fm(tret) 


(/l|/ieob) 


(4.3) 


Erom Eig.m we first reconfirm that the degree of the mis¬ 
match is steeply increased for the last inspiral orbit. This 
indicates that the tidal effect would not be sufficiently 
taken into account in the current version of the EOB 
formalism, although for other inspiral orbits, the perfor¬ 
mance of the EOB formalism appears to be quite good. 
It is also found that the extrapolated waveforms cannot 
be accurately reproduced even if we simply change the 
tidal deformability: If its value is artificially increased, 
the phase evolution is accelerated, and as a result, the 
mismatch is increased in an earlier inspiral stage. If it is 
artihcially decreased, the merger is delayed, and as a re¬ 
sult, the mismatch is badly increased near the last orbit. 


This suggests that a tidal effect, which is not included, 
should be taken into account for improving the perfor¬ 
mance of the EOB formalism. 


V. SUMMARY 

We presented our latest numerical results of longterm 
simulations for the inspiraling binary neutron stars of 
equal mass. By a careful resolution study and extrapo¬ 
lation procedure, we obtain an accurate waveform: The 
estimated total phase error is smaller than ^0.3 radian 
for the total integrated phase of ~ 200 radian and the 
maximum error in the wave amplitude is smaller than 
3%. Using these accurate waveforms, we calibrated the 
waveforms derived by the latest EOB formalism. We 
show that for a binary of compact neutron stars (with 
their radius 11.1km), the waveform by the EOB formal¬ 
ism agrees quite well with the numerical waveform so 
that the total phase error is smaller than 1 radian. By 
contrast, for a binary of less compact neutron stars (with 
their radius 13.6 km), the EOB and numerical waveforms 
disagree with each other in the last a few wave cycles, re¬ 
sulting in the total phase error of ~ 3 radian. We infer 
that this is due to the missing of some tidal effect such as 
nonlinear tidal effect in the current version of the EOB 
formalism, which should be taken into account for im¬ 
proving its performance. 

In this work, we employed only two representative EOS 
and a binary of particular mass. Eor systematically im¬ 
proving the EOB formalism, we have to derive waveforms 
of wider sets of EOS and binary mass. We plan to per¬ 
form more simulations in the future work and to present 
a larger number of the waveforms using the prescription 
developed in this paper. 
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Appendix A: Effective one body formalism 


In this work, we employ an EOB formalism for inspiral¬ 
ing binary neutron stars, which is described in Ref. (2^. 
The base point-particle dynamics for this EOB formalism 
is calibrated by the latest binary-black-hole merger sim¬ 
ulations [i^ and the tidal effects are taken into account 
based on the prescription of Refs. [H, [H, [2l|. Here, we 
briefly review this type of the EOB formalisms and de¬ 
scribe our choice. 

We consider a binary system composed of stars A and 
B with mass of Ma and Mb- The EOB effective metric 
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is defined by 

dslg = -A{r)dt^ + (dd'^ + sin^ 0d(j)^) , 

A{r) ^ ' 

(Al) 


The contribution of tidal effects to the potential is writ¬ 
ten as 

AtidaK?”) = ~ (u) + (A -O' -B)^ , (A12) 

l>2 


where (r, </>) are dimensionless coordinates and their 
canonical momenta are {pr,p^)- We replace the radial 
canonical momentum pr with the canonical momentum 
Pr ,, where a tortoise-like radial coordinate r* is given by 


dr^ _ ^Djr) 
dr A{r) 


(A2) 


Then, the binary dynamics can be described by the EOB 
Hamiltonian 


H^eai{r,Pr,,P<i,) = Mc^-^l + 2v “ l) , (A3) 

where u := MaMb/M'^, M := Ma + Mb, and the effec¬ 
tive Hamiltonian is defined by 


HeS = 


\ 


pI^ + A(r-) |l-|-^-l-2(4-3i^)i^^| . 


(A4) 


The potential A(r) is decomposed into two parts as 


A(r) = App(r) -|- Atidai(?’), (AS) 

where App(r) is the point-particle potential and Atidai(r) 
is the term associated with tidal effects. The point- 
particle potential including up to the fifth PN terms is 

App(r) = pI [l — 2 u + 2uu^ + va^u^ 

+iy{al{i') +a5"lnM)M^ 

-l-(ag(z/) -I- z^ag"(z^) In u)u®] , (A6) 


where u := Ijr, and denotes a (1,5) Fade approx- 
imant. Here, the following coefficients are analytically 
known [13, Ell 


0-4 


al{y) 


a 


In 

5 




94 41 

y “ 

4237 



2 

2275 
512 ” 
41 


64 

T’ 

7004 144 

~to5 


(A7) 

256, 128 

^ (A 8 ) 

(A9) 
(AlO) 


where 7 = 0.5772156 ... is the Euler constant. Following 
Ref. [^, we take the effective form of ag(i^), with which 
results of binary-black-hole-merger simulations are repro¬ 
duced accurately, as 


a^(jy) = 3097.3^ - 1330.6i/ + 81.38. (All) 


where A^^ includes the PN tidal effects and nf is the 
tidal coefficients. Here, the subscripts A and B denote 
the stars A and B. In this work, we include the tidal ef¬ 
fects up to ^ = 4. The coefficient Kf is related to the elec¬ 
tric tidal Love number ki and the compactness C as (see 
Table U for these values of the neutron stars studied in 
this work) 


^ ^ MBMf kf 
' M2'+i ■ 


(A13) 


The tidal potential up to the next-to-next-to-leading cor¬ 
rections is 


The coefficients are 


a 

a 

a 

a 


( 2 ) 

A,1 

( 2 ) 
A,2 

(3) 

AA 

(3) 
A,2 


= 1 + 

analytically known as [l^ 

(A14) 

5a.. 

(A15) 

—A^ + -A^ + 3, 

(A16) 

2 

(A17) 

no 2 311 8 

3 24 3 ’ 

(A18) 


where ■= M^/M. 

Recently, the tidal EOB was improved using resumma¬ 
tion techniques M- We use the gravitational-self-force 
informed I = 2 tidal potential as 


A'Aiu) 


1-f 

-fA 


1 - rLRU 

'^(l-rLR)P’ 


(1 - rLR)V2 


(A19) 


where p is an unknown parameter in the range of 4 < 
p < 6 and we set p to be 4. trr is the light-ring orbit. 
The forms of and are 

= ^J1(1 - aiJi}(l - a2Ji}l-^^^^(A20) 

(A21) 

where the numerical coefficients (ai, 02 , ni, ^ 2 ) are found 
in Ref. M- As in Ref. [l^, we use the tidally corrected 
light ring orbit instead of tlr = 3. For determining the 
value of tlr, we solve the following equation numerically 


A/ ^ 1 dA 

A{u^r)+-ubr — 


= 0 , 


(A22) 


UhR 
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where the tidal part of the potential is included as 
Eq. (IA141) and the value of tlr for the binary neutron 
star models employed in this work is shown in Table I. 
Finally, the potential D[u] v) is given by 


D[u] v) 


1 

1 + Qvu^ + 2(23 — ‘iv)vu'^ 


(A23) 


For calculating the dynamics of the binary orbit under 
the potentials described above, we solve the FOB Hamil¬ 
ton’s equations 


dr 

dt 

d(j) 

dt 

dpr, 

dt 

dptf, 

dt 


A{r) gffreal 

y^D{r) dpr, 

dp^ ’ 

A{r) dH,eal 

^4,- 


(A24) 

(A25) 

(A26) 

(A27) 


Note that we do not include the radial part of the 
radiation-reaction force in Eq. (IA26I) because we 
find this choice advantageous for fitting the extrapolated 
gravitational waveforms. is the radiation-reaction 
force given by 


^ 8 i 

E E (A28) 

1=2 m=l 

where uj = dcji/dt and him denotes the multipolar wave¬ 
forms. Here, him is written as 

him=h^lm + h^m^^^ + hf^^^^^, (A29) 

where h^^ includes the inspiral and plunge waveforms 
given in Ref. [i^, and and are the tidal 

contributions due to the stars A and B. They are given 
by Eqs. (A14)-(A17) of Ref. 
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